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Abstract: An optical analysis is developed to separate forward light scatter 
of the human eye from the conventional wavefront aberrations in a double 
pass optical system. To quantify the separate contributions made by these 
micro- and macro-aberrations, respectively, to the spot image blur in the 
Shark-Hartmann aberrometer, we develop a metric called radial variance 
for spot blur. We prove an additivity property for radial variance that allows 
us to distinguish between spot blurs from macro -aberrations and micro - 
aberrations. When the method is applied to tear break-up in the human eye, 
we find that micro-aberrations in the second pass accounts for about 87% of 
the double pass image blur in the Shack-Hartmann wavefront aberrometer 
under our experimental conditions. 
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1. Introduction 

The anterior corneal surface provides the main optical power for the eye, but its optical 
quality relies on the integrity of the pre-corneal tear film. Normally the tear film varies a few 
microns in thickness across the cornea [1] and provides a smooth optical surface covering the 
microscopically rough corneal surface [2]. The surface topography of the normal tear film, 
therefore, mirrors the macroscopic structure of the corneal surface but hides its microscopic 
structure. Evaporation between blinks causes tear film thinning [3,4] that has the potential to 
change the wavefront aberrations of the eye. Several recent studies using wavefront sensing 
have shown that higher order aberrations (HO A) increase shortly after a blink and, if blinking 
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is suppressed, they can become significantly elevated relative to those observed during the 
normal post-blink period [5-10]. Elevated post-blink increases in HO As have also been 
reported in cases of dry eye [11-14] and in eyes wearing soft contact lenses [15]. These post- 
blink changes in tear film quality are associated with a reduction in optical quality [16-19] 
and loss of visual performance [18,20-22]. Eventually, thinning may cause a localized break 
in the tear film that exposes the irregular, mucin-covered corneal epithelial surface [2], 
introducing small scale optical disturbances that scatter light and further degrade retinal image 
quality [10]. 

Experimental measurements of tear film aberrations have relied primarily upon Shack- 
Hartmann (SH) wavefront aberrometers. However, there is reason to be concerned that SH 
aberrometry lacks the necessary resolution to fully capture the optical disturbances associated 
with the rough corneal surface exposed by tear break-up (TBU) [10,23]. Conventional 
analysis of Shack-Hartmann (SH) images quantifies aberrations based on displacement of a 
multitude of spots formed by an array of lenslets, each of which is the image of a retinal 
beacon of reflected light [24]. Typically the measured spot displacements are subjected to 
modal analysis in which Zernike polynomials are fit to the data by the method of least- 
squares. This modal method of wavefront reconstruction tends to smooth irregularities in the 
data associated with local disturbances in the tear film and may fail to capture most of the loss 
in image quality associated with TBU [23]. Moreover, it is clear from earlier reports of TBU 
[10,23,25], that some SH spots are blurred as well as displaced. These degraded SH spots are 
indicative of poor optical quality within the local region of the eye's optics being sampled by 
an individual lenslet in the aberrometer. Conventional SH wavefront analysis ignores these 
local optical imperfections associated with TBU. 

Our report aims to remedy this shortcoming by introducing a new method of data analysis 
that extracts additional information from the SH spots to quantify wavefront quality locally 
over the domain of each lenslet sampling the eye's pupil. We argue that blurring of SH spots 
is evidence of unresolved micro -aberrations that exist on a fine spatial scale smaller than the 
diameter of the lenslet. This situation is analogous to the inverse scattering problem of 
determining surface characteristics from scattered light measurements [26]. We recover 
information about the magnitude of micro -aberrations from measurements of the size of 
individual SH spots. In this way we extend the utility of the SH aberrometer as a spatially- 
resolved method that also quantifies the micro -aberrations and light scatter at each location in 
the eye's pupil. We note that although TBU provides a motivation for the current study, the 
scope of its applications is not limited to TBU, but is applicable whenever fine micro- 
aberrations are present, for example in cataract [27]. 

Measuring the blur of individual SH spot images is straightforward using a variety of 
image quality metrics [28]. However, in addition to the micro -aberrations present in the region 
of optics represented by each lenslet, these spot images are also blurred by image degradation 
on the first pass through the eye's optics and the resolved macro -aberrations on the second 
pass. Our goal is to apportion these different contributions to spot blur and thereby to quantify 
the magnitude of micro-aberrations present in the eye. To accomplish this goal we develop a 
mathematical model of optical blur produced by a thin random screen using the concept of the 
scattered light field (section 2). To quantify spot blur, we develop a metric called radial 
variance which exhibits a convenient additivity property that enables a simple algorithm to 
estimate the blur produced by each pass in a double -pass SH aberrometer (section 3). We 
validated the algorithm theoretically using local perturbations to the wavefront aberrations and 
demonstrate its feasibility for human eyes (section 4) and investigations of tear film breakup 
(section 5). 

2. Terminology and theoretical background 

In order to describe ocular aberrations with a wavefront error map spanning the full extent of 
the eye's pupil, we sample the emerging wavefront at numerous locations using an array of 
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non-overlapping sub -apertures. For an aberrometer employing a Shack-Hartmann wavefront 
sensor (SHWFS), the subapertures are defined by the faces of individual lenslets in the array. 
If the wavefront over individual subapertures is not flat, local small scale aberrations are 
present which will cause blur in the spot image produced by the lenslet. The classical design 
of the SHWFS assumes the lenslets have sufficiently small diameter that the wavefront is 
approximately planar over each subaperture. In this case, all of the spot images are 
diffraction-limited Airy-disks free of aberration-induced blur and spot displacement is a 
measurement of wavefront tilt over each lenslet. In effect, the wavefront for the whole pupil is 
being approximated by a surface tessellated by flat tiles. The SHWFS determines the slope of 
each tile by measuring the displacement of each spot from the optical axis of the 
corresponding lenslet. The wavefront is then reconstructed from slope measurements 
mathematically by algorithms that either retain all features of the slope data (i.e. zonal 
reconstruction) or that smooth the wavefront to deemphasize local irregularities (i.e. modal 
reconstruction) (Fig. 1). 

Modal reconstructions using Zernike basis functions describe wavefront aberrations by a 
vector of Zernike coefficients associated with a series of Zernike polynomials of increasing 
order. The low-order aberrations (LOA) refer to Zernike polynomials of orders 0-2. The high 
order aberrations (HO A) could, in principle, extend from 3 to N where N is the largest order 
needed to fit the data exactly (i.e. a fully determined least-squares regression). However, in 
practice, HOA are usually specified for a smaller range 3 to R where R < N and frequently R 
< 10. Thus a modal reconstruction of wavefront aberrations based on Zernike coefficients of 
LOA and HOA typically omits the very high order aberrations (VHOA) of orders R to N. The 
total number P of Zernike polynomials contained in orders 0 to R can be determined from the 
recursive sequence P(R) = 1 + R + P(R-l), where P(0) = 1. For example, P(2) = 6. By 
contrast, a zonal reconstruction such as integration by many orbits [29] includes the VHOA. 
Thus the difference between wavefronts reconstructed using modal and zonal methods 
represents the contribution of VHOA, which are known to be indicative of tear film anomalies 
[23,30]. 

order 0 2 10 n 

I I — I — I 
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Macro-aberrations Micro-aberrations 



Fig. 1 . Classification of wavefront aberrations according to the order of Zernike polynomials 
used to represent the wavefront. The magnitudes of the aberration coefficients are not shown 
here. LOA: low order aberrations, HOA: high order aberrations, VHOA: very high order 
aberrations. 

In this paper we use the term "macro-aberrations" to describe wavefronts reconstructed 
from slope measurements. Such aberrations exist on a spatial scale that is large compared to 
lenslet diameter. By comparison, aberrations with a spatial periodicity of oscillation that is 
small compared to lenslet diameter are called "micro-aberrations". Micro-aberrations are 
assumed to be square integrable and may contain random perturbations. They are too fine to 
be resolved by conventional reconstruction methods but nevertheless increase the size of SH 
spots. This division between macro- and micro-aberrations is based on the sampling density of 
the wavefront sensor, analogous to the Nyquist frequency of communication theory that 
separates the frequency spectrum into a band of well-sampled signals and a band of under- 
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sampled signals. In principle, the range of macro -aberrations can be expanded further into the 
domain of micro-aberrations by reducing lenslet diameter. In practice, however, safety 
considerations and increased diffraction by the lenslet aperture limit the density of lenslet 
arrays that can be used with the human eye. For this reason we quantify micro -aberrations by 
taking advantage of the fact that each lenslet in a SHWFS is a miniature optical system that 
produces a blurred image of the retinal beacon. Thus the size of SH spots is a measure of 
wavefront quality over an individual lenslet that is related to the magnitude of micro - 
aberrations. 

To measure micro -aberrations with a SHWFS, we envision fine ripples riding on the 
waves of macro-aberrations (Fig. 2). These ripples (micro-aberrations) add texture to the 
wavefront without affecting curvature or mean slope over the lenslet domain. From a 
geometrical optics perspective, the optical effect of the ripples over any given lenslet is to 
scatter the rays of light in a disorganized plethora of directions, thereby blurring the spot 
image of the retinal beacon [31]. To model such phenomenon mathematically, we include 
multiplicative perturbations in the conventional complex-valued pupil function [32]. This 
extra light field in the traditional pupil function due to micro -aberrations is called the Rytov 
approximation for light scatter in standard references [32]. Since macro -aberrations can also 
blur the SH spot (if they produce sufficient wavefront curvature over individual lenslets), our 
goal was to develop a method for isolating the micro -aberration contribution to blurring of 
individual SH spots. 




Fig. 2. Schematic diagram of an aberrated eye viewing (a) without a random screen, (b) 
through a random screen and sample spot images. Ocular sources of scatter are modeled by a 
thin random screen in the plane of the eye's pupil. The lenslet arrays are colored in orange. 
Courtesy of John. R. Hoffman (Lockheed Martin) during the workshop at Institute for 
Mathematics and Its Applications (IMA) at University of Minnesota. 

Our method of analysis is motivated by the statistical optical description by Goodman [31] 
of imaging through random media. We consider each lenslet of the SHWFS as a small 
imaging system. Phase disturbances over the face of each lenslet are caused by two 
independent sources: gross, structural errors of the eye that cause macro -aberrations and fine, 
micro-aberrations modeled below as a random phase screen. We first briefly review 
Goodman's formulation in the context of the imaging system for human eyes. We then 
develop a metric of radial variance to quantify the blurring of SH spots due to the forward 
propagation of the probe beam that produces the retinal beacon, followed by additional 
blurring due to the reverse propagation of reflected light captured by the SHWFS. 
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2.1 Modeling imaging of the retinal beacon by an individual SHWFS lenslet 

In the imaging optics for the human eye, random phase disturbance can occur from various 
sources. The optical effects of this random disturbance can be statistically modeled, and 
therefore changes of the optical transfer function (OTF) can be modeled by the presence of a 
thin random screen in the optical system [31]. If this thin random screen is located in the pupil 
of the human eye, light fields propagate in a disorganized way and the scattered light fields 
will modify the pupil function P 0 (x,y) by some complex-valued random disturbance t s (x,y). 
We assume the screen is thin enough that the light arrives and leaves the pupil at the same 
location. The new pupil function P(x, y) can be written as the product 

P(x, y) = t s (x 9 y)P 0 (x, y). The optical transfer function for P(x, y) is computed as the 
autocorrelation of the pupil function, 

GO GO 

J J P 0 O, y)P* (x - Xfv x , y - Xfv 2 )t s (x, y)t* (x - Xfv x , y - Xfv 2 )dxdy 
H(v ] ,v 2 ) = =^ . (1) 

V 1 7 l ' GO GO 

J \\P 0 (x,y)\ 2 \t s (x,y)\ 2 dxdy 

X is the wavelength in use, / is the focal length, and (y x , v 2 ) are a spatial frequency pair. The 
asterisk superscript indicates the complex conjugate function. 

Unfortunately, due to lack of knowledge of the transmittance t s (x, y), it is in general 
impossible to analyze Eq. (1) deterministically and so an alternative approach is needed. One 
way to proceed is to assume a certain physical model for t s (x, y) and then derive properties 

of H(v v v 2 ). Alternatively, by assuming that the phase of P 0 follows Gaussian statistics for 
200 eyes, we previously calculated the phase structure function of P 0 and used it to estimate 

the resolution limit (the Fried number) of the eye as an imaging system in terms of the volume 
under the OTF [33]. However, that work used slope measurements from SHWFS and took 
into account only the classical macro -aberrations measured by the SH aberrometer, thereby 
ignoring the micro-aberration contribution over each lenslet. 

Goodman's approach [31] is to look for the ensemble average OTF, i/(v 1? v 2 ), among all 
possible realizations of the random process. The expected value of Eq. (1) is a straightforward 
definition of if(v 1? v 2 ), but since the numerator and the denominator are correlated random 
variables, the computation may not be easy to carry out. Instead, we consider an alternative 
definition for //(v^, v 2 ) as a ratio of expected values, 

g(v„v a )= (2) 
£'[denominator] 

In general, these exact and approximate definitions for H(y l ,v 2 ) are not equal but are 
approximately equal in some cases. For example, if the amplitude of the disturbance function 
is 1, i.e. I t s (x, y) l 2 = 1 for all (x, v), then the denominator in Eq. (1) is a real number and the 
expected values from the two definitions in Eq. (1) and (2) are equal. As a second example, if 
the spatial periodicity of the amplitude I P 0 (x, y) I 2 of the pupil function in the denominator is 

much larger than the correlation width of I t s (jc, y) I 2 then the denominator of Eq. (1) is almost 
constant and independent of each specific realization of the screen. 
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{ { I P 0 (x, y) I 2 1 1, (x, y) I 2 dxdy « J J I P 0 (x, y)\ 2 E[\t, (x, y) I 2 ] dtoiy. (3) 

—00 —00 —00 —00 

Furthermore, if we assume that the transmittance t s (x, y) is a stationary random process, the 
random perturbations are equally likely to occur anywhere in the pupil and the spatial 
correlation functions for the screen t s are independent of x, y. Under this assumption, the 

average OTF for P is written as a product of the OTF, H 0 , for the optical system without the 
screen and the average OTF, H s , for the random screen. 

H(v 1 ,v 2 ) = H 0 (v 1 ,v 2 )H,(v 1 ,v 2 ). (4) 

The equivalent relation for PSF can be obtained by taking the Fourier transform of Eq. (4). 
We use h, h 0 , h s in the lower case for the Fourier transform pairs of H,H 0 ,H s , respectively, 

h = h 0 * h s , (5) 

where h is the average PSF for the whole system (known from the measurement), h 0 is the 
PSF without the random screen (estimated through the standard aberration analysis on the SH 
aberrometer), and h s is the average PSF for the disturbance function due to the random phase 

screen. In summary, the combined blurring effects of micro -aberrations and macro-aberrations 
can be calculated by a simple convolution rule stated by Eq. (5): the combined PSF for any 
given lenslet is the convolution of the PSF due to the micro -aberrations with the PSF due to 
the macro-aberrations. For simplicity, we drop the overline above the average OTF and PSF 
and use the plain letters hereafter. 

Our goal is to obtain information about micro -aberrations from the blurred SH spots in a 
double pass aberrometer in which scatter from micro aberrations occurs when light from a 
narrow laser beam enters the eye, and again within each sub -aperture when light is reflected 
out of the eye. Any additional scattering of light within the retinal tissue that enlarges the 
retinal beacon is bundled with forward scatter in our analysis. Thus the single -pass Eq. (5) is 
employed twice, with possibly imperfect knowledge of the pupil entry point of the laser beam 
on the forward pass. Recovering the full function h s might be possible either using Eq. (4) or 

Eq. (5). However, a method based on Eq. (4) involves dividing H by H 0 could be 
problematic when // 0 has zero values. Also, a method based on Eq. (5) requires a stable 
deconvolution algorithm that is currently lacking. Therefore, rather than attempting the 
difficult problem of recovering the full function h s , we chose instead to recover a scalar 

measure of the width of this spread function. Our measure of choice is the radial variance, 
which has the convenient property of converting convolution into addition, making the 
computation more feasible while seizing a characteristic of the scatter PSF useful for 
characterizing micro-aberrations . 

The model presented above assumes that t s (x,y) is a stationary random process. This 

means that the exact pattern of fine ripples in the wavefront (Fig. 2) are unpredictable, yet 
their statistical properties are stable from one blink to another. Although the magnitude of 
micro-aberrations increases significantly as TBU evolves, and can vary dramatically from one 
part of the cornea to another (see Section 5), it is not unreasonable to assume statistical 
stability at any particular corneal location at some fixed time after blinks occur [34]. 

2.2 Additivity of radial variance 

Although the detailed shape of blurred images can be computed by the convolution 
relationship developed in section 2.1, for the more limited goal of quantifying how much 
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larger an image is compared to the object, it is sufficient to use a scalar measure of size that is 
additive under convolution. This condition is satisfied for univariate density functions by 
variance, defined as the mean-square deviation about the centroid [35]. For bivariate 
distributions of light 7(r,#) in polar coordinates, the corresponding measure of size is radial 
variance V, defined as the second moment of light about the centroid in the radial direction 
r. In Cartesian coordinates, radial variance is computed as 

V(I) = jj^u-uf + (v-vf^I(u,v)dudv, (6) 

where u = jjul(u,v)dudv, v = jjvl(u,v)dudv. A proof of variance additivity for bivariate 

distributions under this definition is included in the Appendix. 

In our optics problem, convolution and variance additivity arise in two contexts. One is the 
standard imaging problem in which the image is the convolution of an object with the PSF of 
the system. Another is the scatter problem in which the average PSF of the whole system is 
computed by convolving the PSF due to macro -aberrations with the average PSF due to a thin 
random screen. Thus we have two variance relations for our imaging systems. First, the radial 
variance of the image is equal to the radial variance of an object plus the radial variance of the 
PSF of the imaging system 

V (Image) = V (Object) + V(PSF^ tem ). (7) 

In the double-pass aberrometer, the object on the first pass is a point source conjugate to the 
retina that produces an image (the retinal beacon) that becomes the object for the second pass. 
The retinal beacon is equal to the PSF of the system from the first pass, including any intra- 
retinal scattering. Second, the radial variance V(h) of the PSF for the eye + lenslet system is 

the sum of the radial variance V(J\) of the PSF for the macro -aberrations over that lenslet 
and the radial variance V(h s ) of the PSF due to the presence of the thin random screen used 
to model micro-aberrations, 

V(h)=V(h,)+V(h s ). (8) 

For the remaining of the paper, we will use these two results repeatedly to estimate the SH 
image blur generated by the eye's optics. 

3. Algorithm for scatter analysis 

In our double-pass optical model, the SH spots are blurred for five reasons: macro and micro 
aberrations encountered on the inward path through a single unknown sub -aperture, retinal 
scatter, and macro and micro aberrations encountered on the second pass through an array of 
lenslets. The macro-aberrations of the eye over each lenslet are known from conventional 
analysis of wavefront slopes but the other factors are unknown. Our goal is to estimate the 
blur due to micro-aberrations on the second pass over each lenslet, despite our incomplete 
knowledge of pupil entry point of the probe beam or the magnitude of retinal scatter. To do 
this we first estimate upper and lower bounds for the radial variance of the retinal beacon. We 
then use those estimates to derive lower and upper bounds for the blur contributed on the 
second pass by micro-aberrations. 

3.1 Upper and lower bounds for the retinal image variance 

In the first pass, the narrow beam of light comes into the eye and forms an image R on the 
retina (Fig. 3). Although the incoming beam does not pass through the lenslet array, we make 
the simplifying assumption that the beam diameter is the same as the diameter of one lenslet. 
In case the beam diameter is different from the lenslet diameter, a telescope system with an 
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appropriate magnification can be used to meet the assumption. A unit-magnification afocal 
telescope makes the lenslet array optically conjugate to the eye's pupil. Thus the pupil entry 
point of the laser beam can be specified in the coordinate system of the lenslet array. We let 
the PSF over the specific subaperture that the probe beam entered be h p . Since the object is a 

point source, V(object) = 0. By using Eq. (7), we conclude that the variance of the retinal 

beacon R is equal to the variance of the PSF for the eye over the sub -aperture defined by the 
entering laser beam plus some additional, unknown amount of variance S contributed by 
retinal scatter, 

V(R)=V(h p ) + S. (9) 

In a carefully aligned measurement system, the position of the lenslet for the entering beam 
relative to the eye's pupil could be known to the experimenter [36]. However, in what follows 
we assume the exact location of the beam entering the eye is unknown. We note that even if 
the pupil entry point of the incoming beam is known, V(R) is unknown due to the lack of 
information about S. 




Fig. 3. A schematic diagram for the double pass optical model for the human eye. R: the retinal 
beacon, h: PSF of the eye + the miniature optical system over each subaperture of the lenslet 
array, h is modeled as a convolution of PSF due to macro -aberrations (h 0 ) and PSF due to 
micro-aberrations (h s ). The SH spot image of R formed by each subaperture is denoted by g. 

In the second pass of light through the eye, the retinal beacon R becomes a common 
object to be imaged by all the lenslets in the aberrometer. Each spot image g.. in the SH 

image is the convolution of the object R with the PSF h.. for the eye + (i, j) -th lenslet that 

quantifies the combined effects of macro- and micro-aberrations. Thus, the variance of each 
spot is the sum of the variances of the retinal beacon and the corresponding eye + lenslet PSF 
(Eqs. (7) and (8)). Specifically, the radial variance of g tj is the sum of the radial variance of 

h.. and the radial variance of R, 

V( gij )=V(h^)+V(R) for each (ij). (10) 

Since V(g tj ) is available from the SH image, solving Eq. (10) for VQiy) requires an estimate 
of V(R). We do this by averaging the lower and upper bounds for V(R) estimated using the 
following rationale. 

Although we cannot measure the retinal beacon directly, the size of the smallest SH spot 
constrains the size of the retinal beacon. For example, since objects are always smaller than 
images, the retinal beacon must be smaller than the smallest SH spot, which means V(R) < 
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min[V(g rc )], where (r,c) is the location of the smallest SH spot. In fact, the retinal beacon must 
be even smaller because some of the variance in SH spots is caused by aberrations and 
diffraction in the exiting path. This blurring in the second pass is minimized when aberrations 
are absent, which means the maximum possible size of the retinal beacon is found by 
assuming the smallest SH spot was produced by a diffraction-limited second pass. This 
argument establishes the upper bound U for V(R) : V(R) < min[V(g rc )] - V(diffraction-limited 
PSF) = : U. 

A lower bound L for V(R) is established by a similar argument. The smallest retinal 
beacon is produced when the probe beam passes through the least aberrated part of the pupil. 
The image of the beacon formed by the lenslet located at the same optimal pupil location will 
produce the smallest SH spot. Thus the smallest SH spot will have been blurred twice by the 
best part of the eye's optical system. Both passes contribute the same amount of radial 
variance to the SH spot, which means the retinal beacon cannot be smaller than half the size 
of the best SH spot. This argument establishes the lower bound LforV(R): V(R) > 
0.5*min[V(g rc )] = : L. 

3.2 Lower and Upper bounds for the PSF blur V(h s ) due to micro-aberrations 

Given these lower bound L and upper bound U for V(R), we can estimate lower and upper 
bounds for V(h) of the PSF created by the second pass through each lenslet. Since the 
following analysis applies to each lenslet, we omit the subscripts for simplicity. Using Eq. 
(10), our estimates for V(h) are as follows: 

V(g)-U<V(h)<V(g)-L, (11) 

where V(g) is the radial variance of the PSF observed in the SH aberrometer. The lower and 

upper bounds for V(h) are represented by two matrices, V(g)-U, andV(g)-L, 

respectively. The amount of the blur due to micro -aberrations can be deduced from these 
matrices as follows. 

We recall Eq. (8) that the radial variance of the PSF h of an optical system is given by a 
sum of the radial variance h 0 of the PSF for macro-aberrations and the radial variance h s of 
the PSF for micro-aberrations: 

V(h)=V(h,)+V(h s ). 

Rearranging the equation, 

VihJ^VW-Vih,). (12) 

Computing the radial variance V(h s ) of the image blur from the micro -aberrations requires 
the radial variance V(h 0 ) of the PSF due to macro-aberrations. The macro-aberrations can be 
estimated by either the modal or the zonal methods and then interpolated over a finer scale to 
create a pupil function for each individual lenslet. We computed the PSF for the macro- 
aberrations using a SH simulation program, details of which are described below in Section 
3.3. 

Finally, in light of Eq. (11), we conclude that the radial variance V(h s ) over each lenslet 
is bounded below by V(g)-U -V(\) and above by V{g)-L-V{h Q )\ 

V(g)-U-V(h 0 )<V(h s )<V(g)-L-V(h 0 ). (13) 
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The lower bound of V(h s ) = V(g)-U —VQiq) is computed with the raw spot image g 
recorded in the aberrometer, the simulated image from macro-aberrations, and the upper 
bound U for V(R)as described above in section 3.1. The upper bound of V(h s ) = 
V(g)-L-V(h 0 ) uses the same information, but with the lower bound L for V(R). We note 
that although V(h s ) varies across different parts of the pupil, the confidence interval for 
V(h s ) is the same as U—L for all lenslets. 

3.3 SH simulation program 

Our SH simulation program was designed to match the specifications of a laboratory SH 
aberrometer used to collect data from human eyes. We use this simulation program to 
calculate the PSF for each portion of the eye's pupil defined by a lenslet, Only macro - 
aberrations are used to calculate the radial variances, V(h DL ) and V(h 0 ). The input for this 

SH simulation program is a high-resolution wavefront map interpolated from measurements 
of wavefront slope. Over each lenslet, we extract the local wavefront phase and use it to 
compute a simulated spot image that has no measurement noise. As an example, in the 
diffraction limited case, our program faithfully represents the Airy pattern in the simulated SH 
image for a point-source retinal beacon. 

Macro-aberrations were represented by wavefronts constructed by either modal or zonal 
methods (Fig. 4). For modal wavefronts, the macro -aberrations are represented as a series of 
Zernike coefficients. The wavefront aberrations can be sampled at any arbitrary rate. The local 
wavefront over each lenslet is unlikely to be planar, and instead contains varying wavefront 
slopes. In the case of the zonal wavefronts, each lenslet produces a single value of wavefront 
slopes. We interpolate these wavefront slopes to generate a planar wavefront surface over 
each lenslet (bi-linear interpolation). A collection of those local planar wavefront surfaces 
provides an approximation for the wavefront map. Therefore, if two neighboring local 
wavefront surfaces happen to have opposite wavefront slopes forming a cusp between them, a 
finite sampling of such wavefronts may lead double spots in one lenslet (Fig. 4). 




(a) (b) 



Fig. 4. An example of the simulated spot images Hq generated by macro-aberrations, (a) the 

modal wavefront, (b) the zonal wavefront. When the wavefront slopes in two neighboring 
lenslets differ and therefore the local wavefront forms a cusp between them. A finite sampling 
of such zonal macro-aberrations may demonstrate double spots in one lenslet (Fig. 4(b)). 

4. Validation using simulated aberrometer images 

In this section, we test the assertion that radial variances of object and PSF can be added 
together to produce the radial variance of the image by using the concept of a localized 
Gaussian disturbance of the wavefront. This validation is placed in the context of tear film 
breakup by modeling the localized thinning of the tear film as an application of a small drop 
of dryness called a "xerop" (from the Greek word xeros (^spoq) for "dry") to the tear layer. 
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The result is a localized shortening of the optical path length from retina to SH wavefront 
sensor that perturbs the wavefront aberration function. For demonstration purposes, we 
assume this perturbation is small enough to fit inside a lenslet face. For simplicity, we choose 
a Gaussian xerop, 



W(x, y) = Cexp 



x 2 + y 2 
2a 2 



(14) 



In the first pass, a narrow beam of light goes through the pupil center when tear film is 
smooth. This first pass optical system is assumed to be diffraction limited in our validation 
test case. This diffraction-limited retinal image from the first pass becomes an object for the 
second pass. Light passing through a xerop on this second pass forms a blurred image in the 
SH image (Fig. 5). 





Fig. 5. Graphical representation of scatter analysis in the double pass setup for Gaussian xerop. 
(a) PSF on the retina from the first pass is the object to be imaged on the second pass. Pixel 
size = 2.72 arcmin. (b) Beam location (1) and several Gaussian xerops (2-4). The Gaussian 
perturbation at location 3 is a drop of wetness that increases optical path length. The Gaussian 
xerops at locations 2 and 4 represent thinning of the tear film that shortens the optical path 
length, (c) PSFs for the second pass. Pixel size = 0.97 arcmin. (d) Double pass SH image. Pixel 
size = 0.97 arcmin. Note that the PSFs in (c) are computed for a point source on the retina. 
Since the retinal image formed from the first pass will contain blur to become an extended 
object for the second pass, the SH images in (d) are not strictly PSFs. For display, a square-root 
transformation was applied to the computed image. 

From these simulated images we computed in Table 1 the radial variances of object (Fig. 
5(a)), PSF (Fig. 5(c)) and images (Fig. 5(d)) associated with Gaussian xerops (Fig. 5(b)). At 
each location, the parameters C and a for Gaussian xerops are also shown in Table 1. For 
all test cases, the error is less than 4% of V(Image). The presence of small errors is due to 
finite discretization of the related functions. 
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Table 1. Parameters and Results of the Virtual Validation using Gaussian Dry Spots 

(Xerops) 



# 


C 


<J 


V(obj) 


V(psf) 


V(obj) + V(psf) 


V(image) 


Error 


1 


1 


0.00 


32.86 


36.57 


69.43 


73.11 


3.68 


2 


1 


33.23 


32.86 


50.12 


82.98 


86.66 


3.68 


3 


-1 


39.88 


32.86 


54.60 


87.46 


91.13 


3.67 


4 


2 


56.97 


32.86 


73.77 


106.63 


110.30 


3.67 



The coefficient C and sigma are in microns. All PSF variances are in arcmin 2 . Xerops are 
sampled at 85x85 points per each lenslet of size 399 micron 2 . The resolution of the SHWFS is 
0.11 arcmin. The wavelength is 633nm. The xerop formed at location 4 is assumed to have a 
negative coefficient, -1, but the radial variance is still positive, error = V(image)-V(obj)- 
V(psf). When sigma = 0, we use a planar wavefront. 

4.1 Resolution of the SHWFS detector 

The additive property of the radial variances is proved analytically in Appendix and is valid 
over the infinite domain. In practice, the additive rule holds only approximately for several 
reasons. The apertures of the optical systems are of finite size whereas the proof in Appendix 
assumes the domain to be of infinite size. The resolution of the SHWFS detector imposes a 
limit since the individual spots are sampled in a discrete grid. Furthermore, the aberrometer 
samples wavefront slope at a finite number of pupil locations. Although each of the factors 
contributes to the total error of the additivity, in this section we focus on detector resolution as 
it applies to the additive property of the radial variances. 

For the simulations shown in Table 1, we assumed that the angular resolution of the 
detector of spot images cast by individual lenslets is 0.1 1 arcmin per pixel. For example, if the 
focal length of the lenslet is 24 mm, this angular resolution is equivalent to a linear resolution 
of 0.76 micron per pixel. Figure 6 shows that the error of the additive rule is about 4% of the 
radial variance V(image) of the image. Since V(image) = V(object) + V(PSF), the error relative 
to V(object) or V(PSF) will be in fact doubled. If the detector resolution decreases, the spot 
image will be more coarsely sampled and the numerical value of V(image) will be less 
accurate. Thus the difference between the numerical and true V(image) will increase, causing 
the error of the additivity to increase. Computer simulations show that as detector resolution 
decreases to 2.7 arcmin per pixel, the error increases to 16.5% of V(image). The error relative 
to V(object) or V(PSF) will be double this value. Furthermore, we observed that the 
computation becomes less reliable. 
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0.18 




Resolution (arcmin) 

Fig. 6. The error of the additive rule of the radial variances with respect to the camera 
resolution. The error was calculated as the difference Viimage) - Vipbject) - V(PSF). As the 
camera resolution decreases, the relative errors of the additive rule increase from less than 2% 
to above 16%. Note that the error was calculated relative to Viimage). Since Viimage) = 
V(object) + V(PSF), the error relative to V(object) or V(PSF) will be 2 times large as the current 
error curve. 

5. Feasibility analysis for human eyes 

5.7 Radial variances of the SHWFS spot images 

The raw SH spot images g were acquired from human eyes by a 633 nm laboratory 
aberrometer [10]. The lenslet diameter was 400 microns referenced in the eye's pupil with the 
telescope magnification of 1. The camera resolution was 6.8 micron per pixel producing a 
59x59 pixel window for each lenslet image. When the tear film forms a smooth surface, 
individual raw spot images are well focused and appear in a regular array as shown in Fig. 
7(a). When the tear film disrupts locally, the spot images become not only displaced from the 
regular pattern, but also blurred due to the light scatter as shown in Fig. 7(b). Tear break-up 
(TBU) was achieved by asking subjects to refrain from blinking for prolonged periods [37]. 
The pupil radius was set at 3 mm for the analysis. Since the lenslet diameter in the SH 
detector is 400 micro meters, the lenslet array contains 15 rows and 15 columns of the spot 
images. The standard aberration analysis from the SH aberrometry with these two arrays of 
the spot images produces two sets of the Zernike aberration coefficients, which are shown in 
Fig. 8. The TBU data exhibits larger higher order aberrations than the base line data, which is 
consistent with the previous literatures [23]. 
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(a) (b) 



Fig. 7. The raw image of the SH spots encircled with a pupil of radius 3 mm. The diameter of 
each lenslet is 400 microns in the SH detector. We expect a 15x15 array of the spot images 
inside the pupil. The ' + ' sign indicates the pupil center, which does not necessarily coincide 
with a spot in a lenslet. The intensity of the raw images was boosted for the display purpose, (a) 
the baseline data when the tear film forms a smooth surface (soon after ablink), (b) the SH 
image after the tear break-up (following blink- suppression). 




(a) (b) 



Fig. 8. The Zernike coefficients for (a) the baseline data and (b) the TBU data. The colors 
represent the WFE. The saturation of the colors represents the magnitude of the Zernike 
aberration coefficients. The TBU data indicate an increase in the amount of HOA. RMS errors 
for the macro-aberrations were calculated with the modal wavefront: (a) 0.50 micron, (b) 1.97 
micron. 

The radial variance maps V(g) of these double pass images are shown in Fig. 9(a) and 

9(b). The original radial variance maps, which hold one value for each lenslet, are bilinearly 
interpolated so that the patterns are more easily readable. Redder colors indicate larger radial 
variances (arcmin 2 ). In general, the spot images possess larger radial variances during the 
TBU (Fig. 9(b)) than for the baseline data (Fig. 9(a)). We argue that since any blur from the 
first pass combined with the retinal scatter will form a common object for the second pass, the 
spatial differences shown in Fig. 9(a) and 9(b) must originate from the second pass. Radial 
variance maps V(h^)of the simulated spot images from second pass macro -aberrations are 

shown in Fig. 9(c) and 9(d) for the baseline data and the TBU data, respectively. We use the 
zonal wavefront to quantify V from macro -aberrations. In contrast to the radial variances of 
the double pass spot images, the radial variances of the simulated spot images for the second 
pass macro-aberrations are in general much smaller in magnitude (median = 13% of the radial 



#139478 - $15.00 USD Received 10 Dec 2010; revised 9 Mar 201 1; accepted 18 Mar 201 1; published 4 Apr 201 1 
(C) 201 1 OS A 1 1 April 201 1 / Vol 19, No. 8 / OPTICS EXPRESS 7431 



variance for the double pass blur for both the baseline and TBU data), which implies that 
macro-aberrations contribute only a small amount to the blur of the raw SH spot images. Their 
overall values and spatial variations across the pupil are small. 




Fig. 9. Comparison of the radial variances V. The units are in arcmin 2 . The original radial 
variance maps contain one value for each lenslet forming maps of 15x15 blocks. We bilinearly 
interpolated the original images so that the patterns are more readily readable. The spot images 
g for the whole eye were acquired by a laboratory aberrometer with the wavelength 633 nm. 
We use the zonal wavefront for the macro-aberrations. The predicted spot images h 0 for the 
resolved aberrations are obtained with the SH simulation program (section 3.3). The radial 
variances for the simulated SH images for the baseline data are close to the DL case, (a) V(g) 
for the base line data, (b) V(g) for the TBU data, (c) V(h 0 ) for the base line data, (d) V(h 0 ) for 
the TBU data. 

5.2 Light scatter for the second pass in a double pass optical system 

The spot blur analysis algorithm in section 3 enables us to isolate the light scattering effect of 
micro-aberrations for the second pass of a double-pass optical system. To demonstrate the 
idea, we apply the algorithm to the data set from section 5.1. A lower bound L and an upper 
bound U for the retinal image quality for the first pass will be computed from image quality of 
the best spot, and therefore include the retinal scatter. The bounds L and U are then averaged 
to calculate the light scatter contributed by the second pass, following Eq. (15). Although any 
of the three criteria - lower bound maps, upper bound maps, or the average maps - for the 
second pass can be equally suitable to quantify the light scatter, we plot the upper and lower 
bound average scatter only for simplicity. The radial variance of the SH spots due to the 
micro-aberrations is isolated by subtracting the blur introduced by the first pass and by the 
macro-aberrations on the second pass from the total blur observed experimentally: 

V(h s )=V(g)-V(hv)-(L+U)/2. (15) 

The values of the radial variance maps for the micro -aberrations are, on average, always 
larger in the TBU data than in the baseline data (Fig. 10). Over the lenslets with large image 
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blur, the average radial variances contributed by the micro -aberrations are considerably 
greater after TBU than before TBU. The radial variances for the micro-aberrations are smaller 
in scale when compared with the radial variances for the double pass total blur as shown in 
Fig. 9(a) and 9(b), but their spatial patterns are very similar because contributions from the 
radial variances V(J\) for the macro-aberrations are almost uniform across the pupil. 

Although the raw spot images for the baseline data (Fig. 7(a)) appear to be well focused, there 
still remains light scatter as indicated in Fig. 10(a). The confidence interval for the scatter can 
be quantified by the difference of the lower and upper bounds for V(h). For each data set, 
although the lower and upper bounds vary across the pupil, as our algorithm in Section 3 
indicates, the confidence interval is uniform. The confidence interval is only slightly larger 
with the TBU data (75.06 arcmin 2 ) than with the base line data (63.3 arcmin 2 ), which indicates 
that the best spots for the two data sets have similar blur in both pre and post TBU. 

The post-TBU radial variance maps in Fig. 10(b) can include the scatter of the eye that 
may not be related to TBU [38]. To see the change in light scatter attributable to TBU, we 
compute the Post-Pre TBU difference maps (Fig. 11). The pattern of blurred spots visible in 
the raw data (Fig, 9b) is similar to the pattern seen in the difference scatter map (Fig. 11). This 
implies that the relative scales of the radial variance maps are independent of the topographic 
characteristics of the patterns. Interestingly, over some regions, the radial variances become 
smaller after TBU (Fig. 11). 




Fig. 10. The average radial variance maps for the light scatter for the second pass only based on 
Eq. (15). (a) V(h . ) for the baseline data. The confidence interval U-L is uniform across the 

pupil and is 63.3 arcmin 2 . (b) V(h ) for the TBU data. The confidence interval U-L is 75.06 

arcmin 2 . The radial variances for the spot images for the TBU data are larger than those for the 
baseline data (e.g. Max = 130.4 arcmin 2 for (a) and 175.3arcmin 2 for (b)). For the baseline data, 
although the raw spot images in Fig. 7(a) appear to be well focused, as Fig. 10(a) indicates, 
there still remains light scatter that has not been resolved by the aberrometer. 
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Fig. 1 1 . The difference of the radial variance maps V (h ) for the light scatter between the 

baseline and TBU data. This map is computed by subtracting Fig. 10(a) from Fig. 10(b) in 
order to see the optical changes of light scatter caused by the event such as TBU. Some regions 
have become negative, which indicates that the light scatter could potentially be reduced during 
the TBU. 

6. Discussion 

6.1 Correlation between the topographic distribution of scatter and the TBU visible in the raw 
SH image 

The spatial distribution of forward light scatter produced by optical perturbations finer than 
can be resolved by our SHWFS (micro -aberrations), was characterized by quantifying spot 
blur in the SHWS images. Since the spot images were acquired in a double-pass optical 
system, the blur due to the first pass and the blur due to the macro -aberrations on the second 
pass have both been subtracted from the measured blur of the raw spot images in order to 
quantify the 2nd pass blur caused by unresolved micro-aberrations (Fig. 10). The result is not 
zero. Although blur from micro -aberrations increases dramatically after TBU, the blur from 
the macro-aberrations on the second pass contributes about 29% of the total blur on the 
second pass for both the baseline and TBU data. 

Similar spatial patterns of the SH spot image blur were reported by Himebaugh et. al [37] 
in cases of TBU. They report that fluorescence and retro -illumination images are indicative of 
the increase in microscopic surface irregularities in regions of TBU. They also show SHWFS 
spots become very blurred in the corresponding TBU regions. By simulation they show that 
the SH spot blur cannot be explained by the second pass macro -aberrations (similar to our 
comparison of Fig. 9(a) and 9(b) to Fig. 9(c) and 9(d)). Therefore, they interpret this as 
evidence that TBU introduces micro-aberrations to the eye's optics. Our study develops a 
metric to quantify the image blur caused by these micro -aberrations. We find that the majority 
of the double pass optical aberrations during TBU are indeed described by the unresolved 
micro-aberrations. From these two studies, we can deduce that the spatial patterns of micro- 
aberrations observed during TBU are likely to be correlated with the spatial distribution of 
exposed rough surface in the tear film. 

To see the changes of the light scatter during TBU, we subtracted the scatter map for the 
baseline data from that for the TBU data (Fig. 11). This difference of the two average scatter 
maps is our preferred measure of the optical change for the following reason: Since the lower 
bound L and the upper bound U for the first pass were estimated based on the quality of the 
best spot in the raw image, both of the bounds include the retinal scatter. The average (L + 
U)/2 is subtracted from the radial variance for the double pass PSF blur. The radial variance 
for the second pass macro-aberrations is further subtracted leaving the radial variance for the 
second pass micro-aberrations. However, it is possible that this scatter map for the second 
pass micro-aberrations might contain the blur from other sources such as the crystalline lens 
or aqueous or vitreous humour, some of which are varying slowly and others of which are 
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varying dynamically during TBU. The difference of the two average scatter maps records the 
optical changes that have occurred dynamically during the intervening time. For a similar 
reason, the difference in the lower bound maps or in the upper bound maps can be used to 
bound such changes. Furthermore, the average maps differ by only a small amount, (U-L)/2, 
from either of the boundary maps, which imply that indeed our bounds are sharp. Comparing 
the difference map in Fig. 1 1 with the raw image in Fig. 9(b) shows that the features in the 
original spot images are also present in the radial variance map. This is a promising optical 
method for understanding time dynamics of micro -aberrations during TBU. The current 
method for spot blur analysis can also be applied to other problems where fine scale 
aberrations are important, e.g. cataract. 

6.2 Resolution of SH aberrometer 

The ability of the analytical method described in this paper to quantify scatter produced by 
micro-aberrations is limited by two instrument-specific parameters. First, the cut-off order for 
the macro -aberrations depends on the resolution of the aberrometer (Fig. 1). Theoretically, as 
we reduce the size of the lenslets, the zonal wavefront may be able to quantify finer scale 
aberrations present in the eye by increasing the cut-off order and therefore expanding the 
range of square integrable functions. 

One currently available high-resolution aberrometer (COAS HD, AMO) has square 
lenslets of size 108 microns and of a focal length 2.421 mm. Its camera captures the spot 
images formed by each lenslet with a 16x16 pixel array. The angular resolution of this system 
for each pixel is therefore 0.108mm/2.421mm x 1/16 x 1807 ;z 'x 60 = 9.58 arcmin/pixel. Our 
analysis from Section 4 suggests that in order to quantify blurring of SH spots, an angular 
resolution of < 1 armin/pixel is required. This is about 1/10 of the angular resolution of the 
high resolution aberrometer that is currently available. The aberrometer used in the current 
study was built in our lab and has a resolution of 0.9740 arcmin/pixel. That is, we have used a 
low resolution aberrometer but with a high sampling density of the SH image, whereas newer 
aberrometers might have higher pupil sampling density, but lower sampling density of the 
image. 

The second important instrument limitation is introduced by the commonly employed 
"range limiting aperture" that prevents highly aberrated rays from being imaged by the 
lenslets. These limiting apertures, located at the common focal point of the relay lenses shown 
in Fig. 3, constrain the angular range of scattered light that will be included within the blurred 
images of the SHWFS spots. Thus largely aberrated rays and scatter with large angles are 
prevented from reaching the lenslets. The current method of spot blur analysis is likely to 
underestimate the scatter of the eye in that case. 

6.3 Micro-aberrations, Strehl ratio and radial variance 

Although the current algorithm quantifies the blur created by the micro -aberrations unresolved 
by the SH aberrometer, it does not quantify the characteristics of these micro -aberrations. 
Some approximate formulas have been suggested by Marechal [39] or Mahajan [40] by 
linking the Strehl ratio to the variance of the total phase aberrations. If the aberrations are 
weak (e.g. RMS of wavefront error < wavelength of light), the overall magnitude of wavefront 
RMS predicts (and can be derived from) the blurred Strehl ratio. 

Alternatively, it is possible to link the Strehl ratio and the micro -aberrations. When we 
decompose the total phase aberrations into the deterministic macro -aberrations and the 
stochastic micro-aberrations, the RMS of the phase aberrations will also be statistically 
varying depending on each realization of the micro -aberrations. The average RMS of the 
phase aberrations can be expressed in terms of the RMS of the macro -aberrations and the 
average RMS of the micro-aberrations. Since the RMS of the macro -aberrations can be 
estimated, we can derive the RMS of the micro -aberrations from the measured Strehl ratio. 
We are currently investigating this approach. 
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Some linkage between the micro-aberrations and the radial variance can also be possible: 
If the phase aberrations are large, the phase difference between two points in the pupil can be 
approximated by the gradient of the phase. Since OTF is the autocorrelation of the pupil 
function, the average OTF can be approximated by the joint characteristic function of the two 
wavefront slopes. The PSF can be represented as the joint probability density function of the 
x- and y- slopes. Since the radial variance of the PSF is simply a sum of two second moments, 
by using the identity in Fourier analysis, it will be equal to the Laplacian of the OTF evaluated 
at the origin. In particular, if we assume that the x- and y- slopes of the wavefronts follow the 
Gaussian distributions with zero means and variances, errand a 2 y , respectively, then the radial 

variance of the PSF will be proportional to the sum cj 2 x +cj 2 . This analysis implies that when 

the phase variance is large, the wavefront slopes undergo large fluctuation. It is the slopes of 
the wavefront that determine the distribution of light energy in the PSF. Applying this 
calculation to the radial variance V(h s ) of the PSF contributed by micro-aberrations, we can 

extract estimates for the variance of the wavefront slopes for the random phase screen, details 
of which will be reported in near future. 

6.4 The whole eye PSF 

Conventional techniques to measure the optical performance of the human eye have been 
dominated by the resolved macro -aberrations. The modal wavefront aberrations are often 
favored. The zonal wavefront aberrations have recently gained popularity among the clinical 
scientists studying tear film optics [23]. However, our results show that across individual 
lenslet apertures the resolved macro-aberrations comprise only 13% of the total aberrations on 
the double pass (or 29% of the total blur created by the second pass) and that much of the blur 
of SH spots, therefore, was caused by the micro -aberrations. This implies that the 
conventional SH analysis is not sufficient to explain all of the optical complications of the eye 
and can be improved by including the micro -aberrations and their blurring of the retinal image 
in order to explain visual performance. Deducing the PSFs of the whole eye from blurring of 
SH spots produced by wave aberrations over individual lenslet apertures is challenging. It is 
unclear whether the PSF for the whole eye can be built from these local PSFs. In order to 
reveal the optical and visual impact of these micro-aberrations, a method for combining these 
spatially uncorrelated (statistical) micro -aberrations from each sub-aperture needs to be 
developed. 

7. Conclusion 

Because of the inherent limitations of SH aberrometer spatial resolution, it is not possible to 
quantify aberrations deterministically beyond the instrument resolution limit. However, we 
have developed a statistical description of the optical contributions of these micro -aberrations. 
Applying this approach to the human eye, we showed that high levels of forward light scatter 
are correlated with the regions of TBU, which we attribute to micro -aberrations. 

Our attention has so far been focused on demonstrating the feasibility of scatter analysis 
algorithm for human eyes. An extended analysis for a large population of eyes or for dry eyes 
could be a useful future study. This method could also be useful to quantify the optical 
qualities of eyes with cataract or refractive surgery or to test the quality of intraocular lenses, 
where fine micro-aberrations are major concerns. Our data confirm that classic SHWFS 
analysis may fail to reveal visually significant optical perturbations, but can be extended 
successfully into the domain of micro -aberrations by the analysis of spot quality. 

Appendix: Radial variance of bivariate random variables and point spread functions 

The mapping between optics and probability theory enables us to explore the radial variances 
both for random variables and for point spread functions (PSF) at the same time. We extend 
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below the optical analysis for the 1 dimensional case by Brace well [35] to the 2 dimensional 
case without imposing additional assumptions on random variables. 

Let X, Y, and Z be an object, an imaging system, and an image, respectively, in two 
dimensional space. 

X = X(u,v), Y = Y(u,v), Z = Z(u,v). 

We denote the point spread function (PSF) for X, 7, Z by p x (u,v), p Y (u,v), p z (u,v), 

respectively. Using standard imaging optics, the PSF for the image is a convolution of the PSF 
for the object and the PSF for the system, 

Pz=Px*Py- ( 16 ) 
A natural choice for generalizing the second moment p z would be 

QO 00 

V(Pz^ = I I [( w_ ^z) 2 +(v -v z ) 2 ~^p z (u,v)dudv, 

-00 -00 

where 

00 00 00 00 

™ z = I I u ^ z v^dudv, v z = I I vp z (u,v)dudv. 

—00 —00 —00 —00 

Since V(p z ) measures the spread of p z in the radial direction, we call V(p z ) the radial 
variance. The definitions for the radial variances for p x and p Y are the same by replacing 
p z with p x and p Y , respectively. We claim that those radial variances add under 
convolution, 



V(p z )=V(Px)+V(p Y \ (17) 



We note that 



V(p z )= J j ( u ~ u z ) 2 p z (u,v)dudv+ j j (v -v z ) 2 p z (u,v)dudv 

—00 —00 —00 —00 

Z , Z —2—2 

= m 20 +m 02 -u z -v z 

m^and m^ 2 are the second order joint moments of the PSF p z (u,v). The similar relation 
holds for p x and p Y . 

Our proof utilizes optical transfer function (OTF). We let O x (o) 1 ,co 2 ), Q Y (a\,a> 2 ), 
O z (a\ , co 2 ) be the OTFs of p x , p Y , and p z , respectively. Since the OTFs and the PSFs are 
Fourier transform pairs, from Eq. (16), we can deduce 

® z (^,^ 2 ) = O x (co [ ,co 2 )O y (co 1 ,co 2 ). (18) 

The OTFs can be expanded in terms of jco x , jco 2 , and moments of the corresponding PSFs. 
The second derivatives of the OTF <& z (co 1 ,a> 2 ) can be found using the second moments 
ofp z : 
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^O z (0,0) = 7 2 m^ ^<J> z (0,0) = /m 0 z 2 . (19) 

OCD x OC0 2 

On the other hand, by applying the product rule to Eq. (18), the second derivatives of 
O z (co 1 ,cd 2 ) are given in terms of the first and second derivatives of O x (co l ,co 2 ) 

and O y (o\ ,co 2 ): 

e 2 ® 7 d 2 o> y _ do Y d® Y _ a 2 o F a 2 o 7 d 2 ® Y _ ^^o Y ao y _ a 2 o y 

z - x -O v +2 — ^ — L + O v f, f = ^O v +2 — ^ — L + O v F 



d&\ 2 d&\ 2 Y dco l dco l x da>l da> 2 da> 2 1 dco 2 dco 2 A dco: 



2 

Evaluating them at (a\ , co 2 ) = (0, 0), 

= f {< + 2« + m 2 y 0 ) , = / (#n£ + 2m 0 > 0 y 1 + m\ 2 ) . (20) 

In fact, Eqns. (19) and (20) are equal. We obtain the following relation, 

m 2 z 0 = m* + 2m^mf 0 + m Y 20 , m 0 z 2 = m* + 2m 0 x 1 m 0 7 1 + m 0 F 2 . 

Similarly, the first moments ra z , m z can be expressed in terms of the first moments for p* v 



and /?f 



Finally, the radial variance for p z 

= (m* + 2mf Q m Y lQ + m Y 20 ) + (m* + 2m*m Y 0l + m 0 7 2 ) - [mf Q + mf 0 )' - (m* + )' 

= ( m 2 X 0 + <2 " (™10 ) 2 " (<1 ) 2 ) + [ m lo + ™02 " {K f ~ {K f ) 

= V(p x ) + V(p 7 ). 
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